#!/usr/bin/perl
use POSIX "fmod"; # float modulus
use List::Util qw[min]; # min function
# GLOBAL
my $grid_size_x = 300;
my $grid_size_y = 300;
my $locitot = 99;
my $name = "sim";
my $vertical_scaling_factor = 1.5; # for hex grid is 0.75, we'll take modulus of 1.5
my $biopsy_size_x = 20; # size of biopsy, 10x10 = 100 cells maximum
my $biopsy_size_y = 10; # biopsies are 1mm by 2mm correspond to 10 crypts by 20 crypts


# Draws random n locations
sub draw_random_locations{
    my ($n,$leftx,$rightx,$topy,$boty) = @_;
    my %h=();
    my $size=0;
    while($size<$n){
        $x = $leftx+int(rand($rightx-$leftx+1)); 
        $y = $topy+int(rand($boty-$topy+1)); 	
	# Convert abs to hex
	my ($hx, $hy) = abs_to_hex_coords($x, $y);
	# Ensure we pick out n DIFFERENT cells
	# do not pick out the same cell twice
        if(!$h{$hx}{$hy}){
            $h{$hx}{$hy}=1;
	    $size++;
        }
    }
    return %h;
}

# convert coords
sub abs_to_hex_coords{
    my ($x, $y) = @_;    
    $y = $y*($vertical_scaling_factor/2.0);
    # Depending on picking y that is divisible by 1.5 or not, add 0.5 to x                                         
    if(fmod($y,$vertical_scaling_factor)){
	$x = $x + 0.5;
    }
    return ($x, $y);
}

sub hex_to_abs_coords{
    my ($x, $y) = @_;    
    # Depending on picking y that is divisible by 1.5 or not, add 0.5 to x                                         
    if(fmod($y,$vertical_scaling_factor)){
	$x = $x - 0.5;
    }
    $y = $y/($vertical_scaling_factor/2.0);
    return ($x, $y);
}

# LOAD pre-computed locations
my $coord, $pn=1;
open (in,"cell.xy.coordinates.per.distance.class.csv");
while(<in>){
    chomp;
    @row=split(/\,/);
    # distance class, pair number, cell 1 x, cell 1 y, cell 2 x, cell 2 y 
    $coord[$row[0]][$pn][0] = $row[1];
    $coord[$row[0]][$pn][1] = $row[2];
    $coord[$row[0]][$pn][2] = $row[3];
    $coord[$row[0]][$pn][3] = $row[4];
    $pn++;
    if($pn==101){$pn=1;}
}
close in;

# Start from sweep 1
# do not use my, so that these hashes are global
 %grid;
 %patterns;
 %indexx, %indexy; 
for $sw ($ARGV[0]..$ARGV[1]){

    # Seed the random number generator with simulation number
    srand ($sw);

    my %grid_coords_x =();
    my %grid_coords_y =();
    open (in,"../$name$sw/$name$sw.crypts");
    while(<in>){
	chomp;
	@row=split(/\,/);
	$grid_coords_x{$row[0]}=$row[1];
	$grid_coords_y{$row[0]}=$row[2];
    }
    close in;

    %patterns=();

    # Load all patterns in memory
    open (in,"../$name$sw/$name$sw.patterns");
    while(<in>){
	chomp;
	@row=split(/\,/);
	$patterns{$row[0]}="";
	for my $loc (0..$locitot){
	    # Hash only the locations that are different from 
	    # the wild type "0" microsattelite length	    
	    # Shaka Zulu proposes hashing only the clone_id key
	    # and as value concatenate "locus,value,locus,value"
	    if($row[$loc+1]!=0){
		# Leading comma ",loc,val,loc,val"
		$patterns{$row[0]}=$patterns{$row[0]}.",".$loc.",".$row[$loc+1];
	    }
	}
    }
    close in;

    
    # Go through the log file of grid state
    %grid=();
    open (i1,"../$name$sw/$name$sw.log");
    while(<i1>){
	chomp;
	@row=split(/\,/);
	# At day 7300, get the grid state
	if($row[2]==7300){
	    # Important - throw it as a double in the hash
	    $grid{$grid_coords_x{$row[0]}*1.0}{$grid_coords_y{$row[0]}*1.0}=$row[1];
	}
    }    
    close i1;    

    open (out3, ">../$name$sw/$name$sw.real.distogram");
    for $dist (1..300){
	for $cell (1..100){
	    $indexx{0}=$coord[$dist][$cell][0];
	    $indexy{0}=$coord[$dist][$cell][1];
	    $indexx{1}=$coord[$dist][$cell][2];
	    $indexy{1}=$coord[$dist][$cell][3];
	    #print STDERR $dist." ".$cell."\n"; #.$indexx{0}." ".$indexy{0}." ".$indexy{1}." ".$indexy{1}."\n";
	    
	    # Compute divergence
	    my ($div,%divpair) = compute_divergence(2,\%grid,\%patterns,\%indexx,\%indexy);
	    my $p = keys(%divpair); # size of hash
	    # If the comparison was between dead cells, do not print it.
	    # so, effectively in each distance class there will be ~98-99 cells, not exactly 100.
	    #print STDERR $p."\n";
	    if($p==1){
		my $genetic_distance = sprintf("%.6f", $divpair{0}{"gd"});
		print out3 $dist.",".$genetic_distance."\n";
	    }
	}
    }
    close out3;
    print STDERR "Sim $sw done\n";
}

    
=pod
# Create allele dropouts
%new_patterns;

# Per panel 60% of 244 loci amplify
my %rem_hash = draw_rand_hash(147,244);
my %index_rem_hash, $ind=0;
for my $key (keys %rem_hash) {
    $index_rem_hash{$key} = $ind;
    $ind++;
}

for my $i (0..($n-1)){
    my $cid = $grid{$indexx{$i}}{$indexy{$i}};
    # Per individual cell, 60% of 147 loci amplify
    my %final_hash = draw_rand_hash(88,147);
    for my $loc (0..243){
	$new_patterns{$cid}{$loc} = $patterns{$cid}{$loc};
	if(!$rem_hash{$loc}){
	    # Globally did not amplify
	    $new_patterns{$cid}{$loc} = -1;
	}
	if(!$final_hash{$index_rem_hash{$loc}}){
	    # Did not amplify in this cell
	    $new_patterns{$cid}{$loc} = -1;
	}
    }
}
=cut


sub compute_divergence{
    # Number of cells
    my ($n,$gridRef,$patternsRef,$indexxRef,$indexyRef) = @_;
    my %grid = %$gridRef;
    my %patterns = %$patternsRef;
    my %indexx = %$indexxRef;
    my %indexy = %$indexyRef;
    my %divergenceofpair = ();

    my $divergence = 0.0;
    my $dead_cells_count = 0;

    my $comparison_count = 0;
    # Compare pattern of cell i to patterns of cell j
    # A total of n*(n-1)/2 comparisons 
    for my $i (0..($n-2)){
	my @col;
	# get the clone_id
	my $cid = $grid{$indexx{$i}}{$indexy{$i}};

	# Skip comparing dead cells
	if($cid==-1){
	    $dead_cells_count++;
	}
	# count if last cell is dead
	if($i==$n-2 && $grid{$indexx{$i+1}}{$indexy{$i+1}}==-1){
	    $dead_cells_count++;
	}

	if($cid!=-1){

	    for my $j (($i+1)..($n-1)){
		my $cid2 = $grid{$indexx{$j}}{$indexy{$j}};

		# Speedup: if cid2 is a dead cell, skip it.
		if($cid2!=-1){
		    # Significant speedup: If the two chosen cells are the same clone
		    # then divergence is 0 between them, thus skip looping over
		    # their microsattelite patterns
		    # Only when they are different clones, compare their pattern
		    my $hamming = 0;
		    if($cid!=$cid2){
						
			# Convert pattern
			my @cell_one = split(/\,/,$patterns{$cid});
			my @cell_two = split(/\,/,$patterns{$cid2});
			
			my %h1 = ();
			my %h2 = ();

			for my $crnch (1..$#cell_one){
			    if($crnch%2==1){
				# Hash locus, value
				$h1{$cell_one[$crnch]} = $cell_one[$crnch+1];
			    }
			}
			for my $crnch (1..$#cell_two){
			    if($crnch%2==1){
				# Hash locus, value
				$h2{$cell_two[$crnch]} = $cell_two[$crnch+1];
			    }
			}

			# Loop over cell one hashed pattern
			while (($locus, $value) = each(%h1)) {
			    #print "$cid $locus $value";
			    # If the locus is present in cell two
			    if(defined($h2{$locus})){
				$hamming += abs($h1{$locus}-$h2{$locus});			    
			    } else {
				$hamming += abs($h1{$locus});
			    }
			}
			#print "Cell TWO:\n";
			# Loop over cell two hashed pattern
			while (($locus, $value) = each(%h2)) {
			    #print "$cid2 $locus $value";
			    # If the locus is present in cell one, do nothing, if absent, compare to 100
			    if(!defined($h1{$locus})){			    
				$hamming += abs($h2{$locus});			 
			    }
			}
			# The rest of the loci are all 100s, so hamming distance is 0.
			# The distance between two cells is hamming/#of valid comparisons
			$divergence += $hamming*1.0/($locitot+1);
			
		    }
		    # Save divergence of the pair of cells whether or not they are from the same clone
		    # Set divergence of the chosen two cells to -1 (comparing dead cells)
		    $divergenceofpair{$comparison_count}{"gd"} = $hamming*1.0/($locitot+1);
		    my ($ax, $ay) = hex_to_abs_coords($indexx{$i}, $indexy{$i});
		    my ($bx, $by) = hex_to_abs_coords($indexx{$j}, $indexy{$j});
		    # Distance base on wrap around hex grid
		    my $xd = min ($grid_size_x-abs($ax-$bx), abs($ax-$bx));
		    my $yd = abs($ay-$by);
		    $divergenceofpair{$comparison_count}{"sd"} = sqrt($xd*$xd+$yd*$yd);
		    $comparison_count++;
		}
	    }
	}
	#print "Comparing cell $i\n";
    }
    #print $n."\n";
    # At the end divide by the total comparisons between cells
    # subtract dead cells count from computation
    # What if all cells that we picked are dead? Then, set divergence to 0.
    if($dead_cells_count<$n-1){
	$divergence = $divergence/(($n-$dead_cells_count)*($n-$dead_cells_count-1)/2.0);
    } else {
	$divergence = 0.0;
    }
    return ($divergence,%divergenceofpair);
}

